
#ifndef AFTERMATH_PROBABILITY_DIST_SCALED_BETA_HPP_INCLUDED
#define AFTERMATH_PROBABILITY_DIST_SCALED_BETA_HPP_INCLUDED

#include "../double_constants.hpp"
#include "dist_traits.hpp"

#include <cmath>
#include <cstdint>
#include <limits>
#include <string>
#include <utility>

namespace aftermath
{
    namespace probability
    {
        // Forward declaration to define traits.
        struct dist_negative_pareto;

        template <>
        struct is_continuous<dist_negative_pareto>
        {
            typedef dist_negative_pareto distribution_type;

            static const bool value = true;
        };

        // A special case of scaled (0; x_max) Beta distribution with parameter beta = 1.
        struct dist_negative_pareto
        {
            typedef dist_negative_pareto type;
            typedef double result_type;

            static const std::string name;

        private:
            double m_alpha, m_x_max;
            double m_cache_expected_value, m_cache_variance, m_cache_standard_deviation;
            double m_cache_axa;

        public:
            dist_negative_pareto()
                : m_alpha(1.0), m_x_max(1.0),
                m_cache_expected_value(0.5),
                m_cache_variance(1.0 / 12),
                m_cache_standard_deviation(double_constants::one_over_root_twelwe),
                m_cache_axa(1.0)
            {
            }

            explicit dist_negative_pareto(double alpha, double x_max)
                : m_alpha(alpha), m_x_max(x_max),
                m_cache_expected_value(x_max * alpha / (alpha + 1.0)),
                m_cache_variance(x_max * x_max / ((alpha + 1.0) * (alpha + 1.0) * (1.0 + 2.0 / alpha))),
                m_cache_standard_deviation(x_max / ((alpha + 1.0) * std::sqrt(1.0 + 2.0 / alpha))),
                m_cache_axa(alpha / std::pow(x_max, alpha))
            {
            }

            double expected_value() const { return this->m_cache_expected_value; }
            double variance() const { return this->m_cache_variance; }
            double standard_deviation() const { return this->m_cache_standard_deviation; }

            double mean() const { return this->expected_value(); }
            double stddev() const { return this->standard_deviation(); }

            double get_alpha() const { return this->m_alpha; }
            double get_x_max() const { return this->m_x_max; }

            double cdf(double x) const { return x <= 0.0 ? 0.0 : (x >= this->m_x_max ? 1.0 : std::pow(x / this->m_x_max, this->m_alpha)); }
            double pdf(double x) const { return (x <= 0.0 || x >= this->m_x_max) ? 0.0 : (this->m_cache_axa * std::pow(x, this->m_alpha - 1.0)); }

            template <std::size_t t_nth_moment>
            double partial_moment(double a, double b) const
            {
                if (b <= 0.0 || a >= this->m_x_max) return 0.0;
                if (b > this->m_x_max) b = this->m_x_max;
                if (a < 0.0) a = 0.0;
            
                // this->m_cache_axa = this->m_xalpha / std::pow(this->m_x_max, this->m_alpha);
                return this->m_cache_axa * (std::pow(b, t_nth_moment + this->m_alpha) - std::pow(a, t_nth_moment + this->m_alpha)) / (t_nth_moment + this->m_alpha);
            }

            bool operator ==(const type& other) const
            {
                return
                    this->m_alpha == other.m_alpha &&
                    this->m_x_max == other.m_x_max;
            }

            bool operator !=(const type& other) const
            {
                return !this->operator==(other);
            }
        };
        
        const std::string dist_negative_pareto::name = "negpareto";
    }
}

namespace std
{
    template <>
    struct hash<aftermath::probability::dist_negative_pareto>
    {
        typedef aftermath::probability::dist_negative_pareto argument_type;
        typedef std::size_t result_type;

        result_type operator()(argument_type const& x) const
        {
            return std::hash<double>()(x.get_alpha()) ^ std::hash<double>()(x.get_x_max());
        }
    };
}

#endif // AFTERMATH_PROBABILITY_DIST_SCALED_BETA_HPP_INCLUDED
